Jamming of hard rods I: From Onsager to Edwards 
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In a 1949 iconic paper, Onsager demonstrated the existence of a first-order isotropic-to-nematic 
transition as the density of an equilibrium suspension of infinitely thin hard rods is increased. When 
the particles in the suspension settle under gravity a new jammed phase emerges. Jammed states 
show reversible behavior when systematically shaken, and therefore, can be treated by statistical 
mechanics as proposed by Edwards. We develop a route to study the phases of jammed packings 
of rodlike particles under the umbrella of Edwards ensemble by extending the framework of Song et 
al., Nature (2008) from hard spheres to hard spherocylinders. In this paper, we start by calculating 
the Voronoi diagram of a set of spherocylinders and present an algorithm for its evaluation. We 
then discuss the next steps to develop the Edwards ensemble in search of novel phases of jammed 
matter of hard rodlike particles. 



I. INTRODUCTION. 

Packings of rods in three dimensions have been exten- 
sively studied due to their rich theoretical [1] and indus- 
trial aspects with applications to lyotropic liquid crystals, 
bio-packings like DNA, powders and construction mate- 
rials 0|- This rich behavior is related to the number of 
degrees of freedom, which gives rise to interesting liquid 
crystal phases beyond the usual crystal-liquid-gas phases. 
However, the existence of intrinsic correlations seems to 
prohibit any theoretical study so far attempted, except 
for infinitely thin hard rods in thermal equilibrium, where 
these correlations vanish as shown by Onsager's virial 
theory Q: Onsager's second order truncation of the free 
energy becomes exact in the limit of infinitely thin par- 
ticles 

When a suspension of rods settles by gravity, a jammed 
phase develops where each particle is at mechanical equi- 
librium and thermal agitation does not play a role. The 
present paper is the starting point of a series that will 
investigate the jammed phases of packings of rods from 
a statistical mechanical point of view using Edwards' en- 
semble theory for jammed granular matter The ap- 
plicability of the theory to the study of same-size ball 
packings in two [5], three a-nd infinite dimensions 

and for polydisperse ball packings [10] motivates the 
present work. As in Onsager's equilibrium theory, rodlike 
jammed particles of different orientations can be treated 
as belonging to different species. Thus, the problem of 
non-spherical monodisperse particles can be mapped to 
the problem of polydisperse spheres and the framework 
developed in [l^l for jamming of different size balls can 
be readily applied to the jamming of hard rods. 

This paper is mainly technical and presents an algo- 
rithm for the calculation of the Voronoi diagram of a set 
of spherocylinders. The Voronoi volume of a single par- 
ticle in Edwards' statistical mechanics (also called the 
volume function [A\) is the analogue of the energy asso- 
ciated to a particle in the classical Boltzmann statisti- 
cal mechanics. We study spherocylinders (a cylinder of 



length L and diameter D = 2a capped with two semi- 
spheres of radius a at both ends) as a model of hard 
rods since their Voronoi diagram is the same as the one 
of line segments. Spherocylinders constitute therefore a 
special case of non-spherical particles which are easier to 
treat than, for instance, ellipsoids of revolution. How- 
ever, the general concepts are expected to apply to these 
more complicated shapes as well. The limiting aspect 
ratio L/D — > oo is discussed as a generalization of On- 
sager's equilibrium results as well as the spherical limit, 
L/D = 0, cLS cl generalization of the equilibrium phases 
of hard spheres. 

The following results are relevant for the present study: 
(i) Onsager's theory for equilibrium thin rods [1, JIj 
L/D — >■ oo, predicts a first order transition between an 
isotropic phase stable below the freezing point at the vol- 
ume fraction 4'\°cc7°'^ L / D ^ 3.29, and a nematic phase 

stable above the melting point 4>\^\t'"^ ^ / ^ ~ 4.19. 

(a) Onsager's transition is analogous to the freez- 
ing transition of equilibrium hard spheres {L/D = 
0) pdl IT^ between an isotropic phase stable below 
(f)'^^'^'^ = 0.494 ± 0.002 and a crystalline phase stable 

above .^^eV"'' = 0-545 ± 0.002. 

(Hi) Philipse et al. experimentally studied isotropic 
jammed phases of thin rods fTsI and showed that the 
volume fraction of the random close packing (RCP) as a 
function of aspect ratio can be fitted as L / D ^ c. 

This formula is valid for L/D ^ 1 and has been related 
to excluded volume arguments from Onsager [131 show- 
ing that c = 5.4 ± 0.2 is half the isostatic coordination 
number of rods, Z{°^ = 10. When L/D->-0 new physics 
beyond the Onsager limit appears as the isostatic limit of 
spherical particles Z^^^ = 6 is approached continuously 
from Zl^^^ M- 

(iv) Jin et al. |l5| found numerically that jammed hard 
spheres undergo a first-order phase transition analogous 
to the equilibrium one. This scenario has been proposed 
by Radin et al. [l6l . \it\ . The jammed isotropic phase 
is stable from random loose packing, RLP, (j\fp~^'^™ ~ 
,sph-jam = ^^Pj^^^^"" ss 0.634 The 



0.536, to RCP, 
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crystal phase is stable from the melting point ''^"^ w 

0.68 to FCC, (/>fP<|'"^™ = tt/VTS « 0.740. Between melt- 
ing and freezing a glass-crystal mixture coexists. 

(v) The application of the Edwards ensemble to 
jammed elongated particles has been treated by Moun- 
fleld and Edwards [18]. This approach uses a phenomeno- 
logical mean-field volume function to describe the mi- 
croscopic interactions between rods. The phase behav- 
ior predicted by Mounfield and Edwards is rather simple 
with no evidence of a discontinuous change in the sym- 
metry of the system, suggesting no phase transition in 
contrast to Onsager's model. We notice that Edwards' 
model, being phenomenological, does not explicitly take 
into account the steric interaction between rods. 

Here we calculate a volume function based on the 
Voronoi volume which accounts for excluded volume ef- 
fects and orientational interactions, crucial to determine 
the phase behavior of the system. The purpose of the 
present work is to investigate the phases of jammed rods 
for all aspect ratios. 




FIG. 1: Voronoi diagram equivalence between two sphero- 
cylinders and two line segments defining the spherocylinders. 
The Voronoi boundary is equivalent to that of two rods of 
vanishing width and length L, i.e., solving for the point where 
Di = D2 is the same as solving Di + a — D2 + a. 



II. EQUIVALENCE BETWEEN THE VORONOI 
DIAGRAM OF SPHEROCYLINDERS AND 
SEGMENTS 

For the study of jammed packings, Edwards' ensemble 
theory uses the volume as a conserved quantity instead 
of the energy as in the classical Boltzmann statistical 
mechanics. The volume associated to each particle is 
the Voronoi volume @ defined as all the points of the 
space that are nearer to the particle than to any other 
one. The Voronoi boundary between two particles is then 
the hypersurface that contains all the points that are 



equidistant to both objects (Fig. [IJ. 

In geometry, the calculation of the Voronoi diagram of 
a set of ellipsoids is an old and complicated problem. The 
calculation of the Voronoi diagram of spherocylinders is 
simpler than the one for ellipsoids thanks to the following 
property: the Voronoi diagram of a set of spherocylin- 
ders of length L and radius a is equivalent to the one 
of segments of length L, see Fig. [TJ This equivalence is 
analogous to the case of equal size spheres. The Voronoi 
diagram of spheres can be calculated by considering only 
the set of sphere centers. Such equivalence disappears for 
instance, when the system consists of spheres of different 
size. 

The radius a of the spherocylinders thus does not ap- 
pear explicitly in the calculation of the Voronoi boundary 
as outlined in the next sections. However, the radius en- 
ters naturally as a limiting condition for the possible con- 
figurations of the spherocylinders: for given orientations 
of the spherocylinders the radius determines a minimal 
separation such that the volumes are not overlapping. 

Furthermore, since the Voronoi diagram of a set of ob- 
jects can be deduced from the one between every pair of 
objects, the problem is reduced to the calculation of the 
Voronoi boundary between two segments. This Voronoi 
boundary is composed of nine analytical surfaces that 
have to be evaluated and properly cut to define the con- 
tinuous and difFerentiable Voronoi boundary. These sur- 
faces arise as the Voronoi boundary due to the following 
interactions which must be treated separately: a point- 
point interaction between any two ends of the segments 
(four interactions in total), a line-point interaction be- 
tween one segment and the end of the other segment 
(four interactions), and a line-line interaction between 
the two segments (one interaction). The principle of the 
algorithm is the calculation of these nine surfaces and 
determine where to cut them to form a continuous and 
differentiable Voronoi boundary. 



III. PARAMETRIZATION OF THE PROBLEM 



We use the parameters displayed in Fig. [5] 
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tively: 



FIG. 2: Parametrization of the problem. 



Let us set the origin of the coordinate system at the 
center of rod i, the vector r points to the center of the 
rod j. We are interested in the calculation of the Voronoi 
boundary along a generic direction s — ss parameterized 
by s. A point on rod i is denoted by the vector ti and a 
point on rod j by tj , where tj originates at the center of 
rod J, i.e., ti and tj denote the directions of rods i and j, 
respectively, such that: ti = tiii, with ti G [—L/2; L/2], 
and analogously for tj . 

In the following we use the convention that a vector 
a can be decomposed as a = ad, where a denotes the 
absolute value and d the unit direction. The product of 
two vectors ab denotes the scalar product a-h — afcfofc, 
where the sum is over all components. 

The Voronoi boundary is obtained from two conditions: 

(i) The point s has the minimal distance to both rod 
i and rod j in the direction s, and 

(a) both distances are the same. 

The square of the distance Di between ti and s and 
Dj between tj and s are 



D^ = {r + t,-s)\ 



(la) 
(lb) 



Condition (i) then requires: 



mi 

dt, 
dti 



0, 



(2a) 
(2b) 



s(Jtj), 



--{s- r)ij = s{si,) - r{fij). 



Condition (ii) requires: 

T-jmin 

i 

which leads to: 

{tf% -s)^^ {t 



7-jmin 
3 ' 
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(3a) 
(3b) 



(4) 



(5) 



The solution of Eqs. © and Q depends on the type 
of interactions between the rods as discussed above. We 
distinguish the four generic cases: 

1. Line-line interaction: i""" £ [—L/2, L/2] and 
if'" G [-L/2,L/2] (one case). 

2. Line-point interaction between the segment i and 
an end-point of j: ij"'" e [-L/2,L/2] and t^ = 
±L/2 (2 cases). 



3. Point-line interaction between the segment j and 
an end-point of i: i™" G [-L/2,L/2\ and t^ = 
±L/2 (2 cases). 



4. Point-point interaction between the end points of i 
and j: ti — ±L/2 and tj — ±L/2 (4 cases). 

By solving Eq. ^ for s, under the conditions that 
s has to be real and positive, one obtains the Voronoi 
boundary at s = ss. The solution of Eq. ([5]) depends 
on the type of interactions between the rods as discussed 
above. We distinguish the four generic cases: 

1. Line-line interaction: i""" G [—L/2, L/2] and 
ij^'" G [-L/2,L/2] (one case). 

2. Line-point interaction between the segment i and 
an end-point of j: i™'" G [-L/2,L/2] and t^ = 
±L/2 (2 cases). 



3. Point-line interaction between the segment j and 
an end-point of i: t™" G [-L/2,L/2] and t, = 
±L/2 (2 cases). 



4. Point-point interaction between the end points of i 
and j: ti = ±L/2 and tj = ±L/2 (4 cases). 

Our algorithm determines the Voronoi boundary be- 
tween two rods as follows: 

1. The configuration of the two rods is specified by 
the directions t^, tj, and the separation r. 

2. Choose a direction s in which the Voronoi boundary 
is to be determined at s — ss. The value of the 
boundary s will be a function of the configuration 
of the rods and s: 



This leads to the minimal values along ti and tj , respec- 



s{ti ,ij,S,T^. 



(6) 
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3. Solve Eq. (O for each of the nine interaction cases. 
If the resulting s is real and positive it is a valid 
solution for the Voronoi boundary in the direction 
s. 

It is important to note that depending on the separa- 
tion f and orientation of the two rods it is possible that 
two, or more, different types of interaction contribute to 
the Voronoi boundary for a given direction s. This is 
illustrated in Fig. [2l where in the direction s' both line- 
line interactions and point-point interactions contribute. 
This ambiguity is considered in the algorithm, because 
for multiple contributing interactions one simply obtains 
multiple valid solutions for s. 



IV. SOLUTION FOR THE DIFFERENT 
INTERACTIONS 

We start by choosing a given direction s along which 
we will calculate the Voronoi boundary specified by s. 
The relevant scalar products between the unit vectors, 
s, f, ti and tj can be obtained. We consider the first 
interaction between the rods. 



A. Line-line interaction 

This case arises if the two solutions ([3]) fall inside the 
length of the segments. The conditions are: 



e[~L/2,L/2]. 



j.inin 



(7a) 
(7b) 



In this case t^*" and t™*" are given by Eqs. ([3]). Sub- 
stituting these expressions into Eq. ^ then leads to a 
quadratic equation for the value s of the boundary: 

.2 r 1 s r 1 

{sU)'^-{sijf +2- {sij){fij)~fs +l-{rljf ^Q. 

(8) 

Setting I = s/r, the Voronoi boundary between the 
two rods thus scales with the separation r in this case: 



rl{ii,ij,f,s), 



(9) 



where I is obtained from Eq. ([8]) and only depends on 
the four different directions. 



B. Line-point interaction 

In this case the solution for i™'" falls along the line 
segment i and t™'" is at one of the end points of rod j. 

We choose the top of tj as the point (the other case is 
analogous) and we obtain: 



if" ^s{sU) e [-L/2,L/2] 
^ " 2 ' 



(10a) 
(10b) 



the second equation arises since 

s(%) - r(ff,) ^ hL/2,i/2]. (11) 

Substituting the expressions Eq. (fTO|) into Eq. ([5]) then 
leads to a quadratic equation for s: 



(12) 



C. Point-line interaction 

This interaction is analogous to line-point. We show 
the case when the point is at the top of U. The conditions 
are: 



2' 



tf'^^s{si,)~riH,)e[-L/2,L/2]. 
The first equation arises since 



(13a) 
(13b) 



s{sU) (^[-L/2,L/2]. (14) 
The quadratic equation resulting from Eq. ([S]) is then: 



sL,^., (L 

{sU)+[ — 

r r \2r 



(fs) - {sij){fij) 



0. (15) 



D. Point-point interaction 

In this case the conditions are: 

s{sU)i[-Ll2,Ll2l 
s{sij) ^ r{H^) i [-i/2,L/2], 



(16a) 
(16b) 



so that the two points if" and i™'" are both fixed and 
equal to L/2 or —L/2. 
We choose the values 



C" = L/2, 
if" = L/2. 

Solving Eq. © with ((T71), we find: 



2{fs) + H{st,)-{st,)) 



(17a) 
(17b) 



(18) 



The two line-point interactions and the point-point inter- 
action do not scale with r except in the limit L/r 0, 
in this case from the point-point solution Eq. (|18p we re- 
cover the Voronoi boundary between two hard spheres as 
calculated in (1, 01 : 

(19) 



2{rs)' 



5 



EXAMPLES OF SOLUTIONS AND 
DISCUSSION 



but finite rods, Cu will intervene much more than point 
interactions. Thus, 



As an example of the algorithm, we apply it to dif- 
ferent situations in 2 dimensions. Fig. [3] We consider 
the rod i on the left and the rod j on the right. The 
top-left panel shows the Voronoi boundary in different 
colors corresponding to different interactions: a point- 
point interaction at the top of the boundary in red, then 
a point-line interaction in green, then a line-line interac- 
tion in blue, then another point-line interaction in green 
and so on. The other panels are analogous. 



Cpp < C'lp < C;/, 



(20) 



and reverse for small rods. We expect that: 

(i) For very large ratio only Ca needs to be considered. 

(ii) For large ratio, only Cu and Cp/. 
(Hi) For small ratio, only Cpp and Cpi. 

(iv) For very small ratio (slightly deformed spheres), 
only Cpp. 

These remarks are important to understand the be- 
haviour of the volume fraction as a function of aspect 
ratio and will be analyzed in future works. 





FIG. 3: Voronoi boundary in 2d for difTerent situations. Each 
color of the boundary represents a type of interaction. Red = 
point-point, green = point-line and blue = line-line. 

An important fact we should underline is that the 
three types of interactions are each composed by four 
(for point-point and point-line) or one (for line-line) of 
the nine surfaces. In dimension greater than 2 and for 
general angular configurations, the equation describing 
the surfaces of the point-line interaction and the line-line 
interaction are quadratic and the ones of the point-point 
interaction are linear. We thus obtain a surface made by 
a selection of four planes (due to the point-point inter- 
action), four paraboloids of revolution (point-line) and 
one hyperbolic paraboloid (line-line) cut and paste to- 
gether. In the limit L — > oo, we only keep the hyperbolic 
paraboloid. The two dimensional case is special because 
the line-line equation always gives birth to a line. We 
thus obtain a curve made by a selection of five lines and 
four parabola. 

The contribution of the different interactions, noted: 
Cu for line- line, Cpi for point-line and Cpp for point-point, 
naturally changes with the ratio L/D and may lead to 
different trends in the behaviour of the volume fraction 
of the packing as a function of the aspect ratio, thus 
explaining the results of It is evident that for long 



VI. CONCLUSION 

We have shown how to calculate the Voronoi boundary 
between two spherocylinders as a boundary composed by 
analytical surfaces. The calculation of the Voronoi dia- 
gram for a set of spherocylinders is then possible as well 
as the calculation of the Voronoi volume. The analytical 
expression for the Voronoi boundary is highly simplified 
in the limit L — >■ oo, i.e., in the Onsager limit for finite 
radius a, since only the solution scaling with r given by 
the line-line interaction Eq. ([5]) needs to be considered; 
the surface is then a hyperbolic paraboloid. The differ- 
ent solutions involving point interactions appear because 
of the finite extension of the spherocylinders and make 
the boundary more complicated and not scaling with r. 
These cases do not appear for infinite extension of the 
rods, thus Onsager limit corresponds to the simplest one. 

The next step is to incorporate these expressions in or- 
der to calculate the excluded volume and excluded sur- 
face for two interacting particles. The statistical theory 
for the mean Voronoi volume of a rod should be similar 
to that of (multi disperse) spheres (T^. Thus, the cal- 
culation of the probability distribution to find a particle 
contributing to the Voronoi boundary can be performed. 
This calculation is analogous to the sphere case and can 
be performed by integrating over Theta functions as in 
[Q-Q . In addition to the integration over the vector r one 
eventually also has to integrate over the possible direc- 
tions of the rod j, tj. Such a calculation is in progress. 

It is interesting to compare the case of rods of finite 
length to the one of infinite length. The thickness D ~ 2a 
of the rods does not intervene in any of the formulas of 
the Voronoi boundary in both cases and will intervene 
only to determine the possible interactions. For spheres, 
we have approximated [9.] the pair correlation function 
to a step function for the background term and a delta 
function for the contact term and, on average, this ap- 
proximation has small influence on the final prediction of 
the volume fraction of the packing. For rods with infi- 
nite length, the pair correlation function (as a function 
of distance and angles) is equiprobable for non overlap- 
ping configurations as shown in [3| . This means that the 
approximation of ^, i9|] for the background and contact 
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terms may becomes exact for infinite rods. For rods of ble, and considering the pair correlation function as a 
finite length, the correlation function is not equiproba- step-function becomes an approximation. 



G. J. Vroege and H. N. W. Lekkerkerker, Rep. Prog. 
Phys. 55, 1241 (1992). 

P.-G. de Gennes. The physics of liquid crystals (Oxford, 
Clarendon, 1974). 

L. Onsager, Ann. N. Y. Acad. Sci. 51, 627 (1949). 

S. F. Edwards and R. B. S. Oakeshott, Physica A 157, 

1080 (1989). 

S. Meyer, C. Song, Y. Jin, K. Wang and H. A. Makse, 
Physica A 389, 5137 (2010). 

C. Song, P. Wang and H. A. Makse, Nature 453, 629 
(2008). 

C. Song, P. Wang, Y. Jin, H. A. Makse, Physica A 389, 

4497 (2010). 

P. Wang, C. Song, Y. Jin, H. A. Makse, Physica A 390, 
427 (2010). 

Y. Jin, P. Charbonneau, S. Meyer, C. Song and F. Zam- 
poni, Phys. Rev. E 82, 051126 (2010). 
M. Danisch, Y. Jin and H. A. Makse, Phys. Rev. E 81, 
051303 (2010). 

W. W. Wood and J. D. Jacobson, J. Chem. Phys. 27, 
1207 (1957). 



[12] B. J. Alder and T. A. Wainwright, J. Chem. Phys. 27, 

1208 (1957). 
[13] A. P. Philipse, Langmuir 12, 1127 (1996). 
[14] A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Still- 

inger, R. Connelly, S. Torquato, P. M. Chaikin, Science 

303, 990 (2004). 
[15] Y. Jin and H. A. Makse, Physica A 389, 5362 (2010). 
[16] C. Radin, J. Stat. Phys. 131, 567 (2008). 
117] D. Aristoff and C. Radin, J. Math. Phys. 51, 113302 

(2010). 

]18] C. C. Mounfield and S. F. Edwards, Physica A 210, 279 

(1994). 

119] J. Blouwolff and S. Fraden, Europhys. Lett. 76, 1095 
(2006). 



Acknowledgements. This work is supported by the 

National Science Foundation and Department of Energy. 
The authors would like to thank Jean-Marc Mercier and 
Yuliang Jin for helpful discussions. 



